%————————————————2006。8。9——————————————————
%———————————————creat by liu——————————————————
%——————————create for 8CPM+TCM text by myself———————————
%——————————————TCM parameters——————————————————
%—————编码多项式：Polynomial Description of code (15,2)—————————
%———The length of a minimum distance event is 3.2822_欧式距离?——————
%——————————The minimum decision depth is ?—————————————
diary TCM152_M8_h8_1RC_quasi_gozero3.txt
symbol_number=46;
id=sqrt(-1);
pi=3.14159;
h=1/8;                                       % 调制指数  
M=8;                                         % 进制数  
sample_number=8;                            %每码元样点个数
Tb=1/32000;                                  % 码元持续时间
fs=1/Tb*sample_number;                       % 基带取样率
tt=0;
BER_SNR=zeros(10,1);
for SNR=7:0.5:7
    tt=tt+1;
    it_num=5000;
    BER=zeros(1,it_num+1);
    BER(1)=0;
    for it=2:it_num+1


        % 
        rate=2/3;
        cl=4;                                            
        gen=[15 2];                                       
        memory=cl-1;                                     
        trellis_cc=poly2trellis(cl,gen);   

        map=[-7,-3,1,5,
             -3,-7,5,1,
             -5,-1,3,7,
             -1,-5,7,3,
             -3,-7,5,1,
             -7,-3,1,5,
             -1,-5,7,3,
             -5,-1,3,7];
            
         state_trans_F=[1,5,1,5,
                        1,5,1,5,
                        2,6,2,6,
                        2,6,2,6,
                        3,7,3,7,
                        3,7,3,7,
                        4,8,4,8,
                        4,8,4,8];                
                                            
        bit_msg=randsrc(1,symbol_number*2,[0,1]);
        bit_msg(1:6)=[0 0 0 0 0 0]; bit_msg(44*2-1:46*2)=[0 0 0 0 0 0];
        for t=1:symbol_number
            msg(t)=2*bit_msg(2*t-1)+bit_msg(2*t);
        end
        state_temp=1;
        for t=1:symbol_number
            code(t)=map(state_temp,msg(t)+1);
            state_temp=state_trans_F(state_temp,msg(t)+1);
        end
        Ik=code;
        %-----------------------------------------------------------

        %-----------------------------------------------------------
       
        L=1;
        g=zeros(1,sample_number+1);
        for i=0:L*sample_number
            t=i/fs;
            g(i+1)=1/(2*L*Tb)*(1-cos(2*pi*t/(L*Tb)));
        end;

        q=zeros(1,L*sample_number+1);
        for i=1:L*sample_number+1
            t=i/fs;
            q(i)=t/(2*L*Tb)-1/(4*pi)*sin(2*pi*t/(L*Tb)); 
            %q(i)=t/(2*L*Tb); 
        end;

        %---------------------------------------------------------

        %--------------------------------------------------------
        % 计算载波相位  
        % 模2*pi的绝对相位值（弧度）
        theta=zeros(1,symbol_number*sample_number);                               
        I_signal=zeros(1,symbol_number*sample_number); 
        Q_signal=zeros(1,symbol_number*sample_number);

        carrier_phase=pi/2;

        sigma_a=zeros(1,symbol_number);
        fi=zeros(1,symbol_number*sample_number);
        fi1=zeros(1,symbol_number*sample_number);
        for n=1:symbol_number
            if n==1
               sigma_a(n)=0;  
            else
               sigma_a(n)=sigma_a(n-1)+Ik(n-L); 
            end
            theta_n=sigma_a(n)*h*pi;
            delta=zeros(1,sample_number);
            for m=1:sample_number
                delta(m)=2*pi*h*Ik(n)*q(m);
                fi(sample_number*(n-1)+m)=delta(m)+theta_n;
                fi1(sample_number*(n-1)+m)=delta(m)+theta_n+carrier_phase;
            end
        end
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

        I_signal=cos(fi);              % 基带I路信号
        Q_signal=sin(fi);              % 基带Q路信号    
        s=I_signal+Q_signal*id;

        I_signal1=cos(fi1);            % 基带I路信号 有相偏
        Q_signal1=sin(fi1);            % 基带Q路信号 有相偏   
        s1=I_signal1+Q_signal1*id;

        N=130;
        F=[0,27500,29000,fs/2]*2/fs;
        A=[1,1,0,0];
        lpf=firls(N,F,A);
        I_s_lpf=conv(I_signal1,lpf);
        Q_s_lpf=conv(Q_signal1,lpf);
        I_signal11=I_s_lpf(N/2+1:length(I_s_lpf)-N/2);
        Q_signal11=Q_s_lpf(N/2+1:length(I_s_lpf)-N/2);
        
        sigma=sqrt(sample_number/(2*log2(M)*10.^(SNR/10)));
        I_noise=zeros(1,symbol_number*sample_number);
        Q_noise=zeros(1,symbol_number*sample_number);
        for i=1:2:symbol_number*sample_number
              [I_noise(i) I_noise(i+1) ]=gngauss(sigma);
              [Q_noise(i) Q_noise(i+1) ]=gngauss(sigma);
        end;
        
        N=130;
        F=[0,0.6/Tb,0.7/Tb,fs/2]*2/fs;
        A=[1,1,0,0];
        lpf1=firls(N,F,A);
        I_lpf=conv(I_noise,lpf1);
        Q_lpf=conv(Q_noise,lpf1);
        I_noise=I_lpf(N/2+1:length(I_lpf)-N/2);
        Q_noise=Q_lpf(N/2+1:length(I_lpf)-N/2);
        r=I_signal11+I_noise+(Q_signal11+Q_noise)*id;
        %r=I_signal1+I_noise+(Q_signal1+Q_noise)*id;

        Z1=sum(r(1:8)*conj(s(1:8)).');
        fi_est=atan2(imag(Z1),real(Z1));

        %---------TCM decode-------------------------------------------------------
        % viteri 检测
        alpha=0.005;
        symbol_number=length(r)/sample_number;
        state_number=8;
        len=length(r)/sample_number;
        state_trans_B=[1,1,2,2,
                       3,3,4,4,
                       5,5,6,6,
                       7,7,8,8
                       1,1,2,2,
                       3,3,4,4,
                       5,5,6,6,
                       7,7,8,8];
        Un_trans=[-7,1,-3,5,
                  -5,3,-1,7,
                  -3,5,-7,1,
                  -1,7,-5,3,
                  -3,5,-7,1,
                  -1,7,-5,3,
                  -7,1,-3,5,
                  -5,3,-1,7];

        initial_metric=[0,-999*ones(1,state_number-1)];
        metric_state=-99999999*ones(state_number,symbol_number);
        sum_Uk_old=zeros(state_number,1);
        sum_Uk_new=zeros(state_number,1);
        msg_decision=zeros(state_number,symbol_number);               
        temp_msg_decision=zeros(state_number,symbol_number);  
        last_Vn=zeros(1,state_number);

        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        for n=1:len
          r_T=r((n-1)*sample_number+1:n*sample_number)*(cos(fi_est)-id*sin(fi_est));
           %*******************************************************************************************************************
          if n==1
            %@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ 
            for i=1:state_number
                for j=1:4
                    Vn=sum_Uk_old(state_trans_B(i,j));
                    Un=Un_trans(i,j);
                    refe_fi=zeros(1,sample_number);
                    for k=1:sample_number
                      refe_fi(k)=pi*h*Vn+2*pi*h*Un*q(k);
                    end 
                    ss=cos(refe_fi)+id*sin(refe_fi);
                    metric=real(sum(r_T*conj(ss).'));
                    if metric+initial_metric(state_trans_B(i,j))>metric_state(i,n)
                        metric_state(i,n)=metric+initial_metric(state_trans_B(i,j));
                        Uk_decision(i,n)=Un_trans(i,j);
                        sum_Uk_new(i)=Vn+Un;
                        last_Vn(i)=Vn;
                        msg_decision(i,:)=temp_msg_decision(state_trans_B(i,j),:);
                        msg_decision(i,n)=fix((i-1)/4)+2*rem((j-1),2);
                    end
                end
            end
            temp_msg_decision=msg_decision;
  
            ss=s((n-1)*sample_number+1:n*sample_number);  
            en_1=imag(sum(r_T.*conj(ss))); 
            fi_est=fi_est+alpha*en_1;                          % 更新载波相位估计值
            %@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

          elseif n<=43 
             for i=1:state_number
                for j=1:4
                    Vn=sum_Uk_old(state_trans_B(i,j));
                    Un=Un_trans(i,j);
                    refe_fi=zeros(1,sample_number);
                    for k=1:sample_number
                      refe_fi(k)=pi*h*Vn+2*pi*h*Un*q(k);
                    end 
                    ss=cos(refe_fi)+id*sin(refe_fi);
                    metric=real(sum(r_T*conj(ss).'));
                    if metric+metric_state(state_trans_B(i,j),n-1)>metric_state(i,n)
                        metric_state(i,n)=metric+metric_state(state_trans_B(i,j),n-1);
                        Uk_decision(i,n)=Un_trans(i,j);
                        sum_Uk_new(i)=Vn+Un;
                        last_Vn(i)=Vn;
                        msg_decision(i,:)=temp_msg_decision(state_trans_B(i,j),:);
                        msg_decision(i,n)=fix((i-1)/4)+2*rem((j-1),2);
                    end
                end
             end
             temp_msg_decision=msg_decision;
            %@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
            
            if n<=3
             ss=s((n-1)*sample_number+1:n*sample_number);   
             en_1=imag(sum(r_T.*conj(ss))); 
             fi_est=fi_est+alpha*en_1;                          % 更新载波相位估计值 
            else
             max=1;
             for i=2:state_number
               if metric_state(i,n)>metric_state(max,n)
                     max=i;
               end
             end
             last_phase_state=last_Vn(max);
             Uk_sur=Uk_decision(max,n);
             refe_fi=zeros(1,sample_number);
             for k=1:sample_number
                refe_fi(k)=pi*h*last_phase_state+2*pi*h*Uk_sur*q(k);
             end 
             ss=cos(refe_fi)+id*sin(refe_fi); 
             ss1=s((n-1)*sample_number+1:n*sample_number);   
             en_1=imag(sum(r_T.*conj(ss))); 
             fi_est=fi_est+alpha*en_1;                          % 更新载波相位估计值 
            end


         %@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

          elseif n==44 
             for i=1:4
                for j=1:2:3
                    Vn=sum_Uk_old(state_trans_B(i,j));
                    Un=Un_trans(i,j);
                    refe_fi=zeros(1,sample_number);
                    for k=1:sample_number
                      refe_fi(k)=pi*h*Vn+2*pi*h*Un*q(k);
                    end 
                    ss=cos(refe_fi)+id*sin(refe_fi);
                    metric=real(sum(r_T*conj(ss).'));
                    if metric+metric_state(state_trans_B(i,j),n-1)>metric_state(i,n)
                        metric_state(i,n)=metric+metric_state(state_trans_B(i,j),n-1);
                        Uk_decision(i,n)=Un_trans(i,j);
                        sum_Uk_new(i)=Vn+Un;
                        last_Vn(i)=Vn;
                        msg_decision(i,:)=temp_msg_decision(state_trans_B(i,j),:);
                        msg_decision(i,n)=fix((i-1)/4)+2*rem((j-1),2);
                    end
                end
             end
             temp_msg_decision=msg_decision;
            max=1;
            for i=2:state_number
               if metric_state(i,n)>metric_state(max,n)
                     max=i;
               end
            end
            last_phase_state=last_Vn(max);
            Uk_sur=Uk_decision(max,n);
            refe_fi=zeros(1,sample_number);
            for k=1:sample_number
                refe_fi(k)=pi*h*last_phase_state+2*pi*h*Uk_sur*q(k);
            end 
            ss=cos(refe_fi)+id*sin(refe_fi);  
            ss1=s((n-1)*sample_number+1:n*sample_number);   
            en_1=imag(sum(r_T.*conj(ss))); 
            fi_est=fi_est+alpha*en_1;                          % 更新载波相位估计值 
            %@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
            
            %@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

           elseif n==45 
             for i=1:2
                for j=1:2:3
                    Vn=sum_Uk_old(state_trans_B(i,j));
                    Un=Un_trans(i,j);
                    refe_fi=zeros(1,sample_number);
                    for k=1:sample_number
                      refe_fi(k)=pi*h*Vn+2*pi*h*Un*q(k);
                    end 
                    ss=cos(refe_fi)+id*sin(refe_fi);
                    metric=real(sum(r_T*conj(ss).'));
                    if metric+metric_state(state_trans_B(i,j),n-1)>metric_state(i,n)
                        metric_state(i,n)=metric+metric_state(state_trans_B(i,j),n-1);
                        Uk_decision(i,n)=Un_trans(i,j);
                        sum_Uk_new(i)=Vn+Un;
                        last_Vn(i)=Vn;
                        msg_decision(i,:)=temp_msg_decision(state_trans_B(i,j),:);
                        msg_decision(i,n)=fix((i-1)/4)+2*rem((j-1),2);
                    end
                end
             end
             temp_msg_decision=msg_decision;
             
            max=1;
            for i=2:state_number
               if metric_state(i,n)>metric_state(max,n)
                     max=i;
               end
            end
            last_phase_state=last_Vn(max);
            Uk_sur=Uk_decision(max,n);
            refe_fi=zeros(1,sample_number);
            for k=1:sample_number
                refe_fi(k)=pi*h*last_phase_state+2*pi*h*Uk_sur*q(k);
            end 
            ss=cos(refe_fi)+id*sin(refe_fi);  
            ss1=s((n-1)*sample_number+1:n*sample_number);   
            en_1=imag(sum(r_T.*conj(ss))); 
            fi_est=fi_est+alpha*en_1;                          % 更新载波相位估计值 
            %@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

            elseif n==46 
             for i=1:1
                for j=1:2:3
                    Vn=sum_Uk_old(state_trans_B(i,j));
                    Un=Un_trans(i,j);
                    refe_fi=zeros(1,sample_number);
                    for k=1:sample_number
                      %refe_fi(k)=2*pi*h*Vn+4*pi*h*Un*q(k)+W(k);
                      refe_fi(k)=pi*h*Vn+2*pi*h*Un*q(k);
                    end 
                    ss=cos(refe_fi)+id*sin(refe_fi);
                    metric=real(sum(r_T*conj(ss).'));
                    if metric+metric_state(state_trans_B(i,j),n-1)>metric_state(i,n)
                        metric_state(i,n)=metric+metric_state(state_trans_B(i,j),n-1);
                        Uk_decision(i,n)=Un_trans(i,j);
                        sum_Uk_new(i)=Vn+Un;
                        last_Vn(i)=Vn;
                        msg_decision(i,:)=temp_msg_decision(state_trans_B(i,j),:);
                        msg_decision(i,n)=fix((i-1)/4)+2*rem((j-1),2);
                    end
                end
             end
             temp_msg_decision=msg_decision;
             
            max=1;
            for i=2:state_number
               if metric_state(i,n)>metric_state(max,n)
                     max=i;
               end
            end
            last_phase_state=last_Vn(max);
            Uk_sur=Uk_decision(max,n);
            refe_fi=zeros(1,sample_number);
            for k=1:sample_number
                refe_fi(k)=pi*h*last_phase_state+2*pi*h*Uk_sur*q(k);
            end 
            ss=cos(refe_fi)+id*sin(refe_fi);  
            ss1=s((n-1)*sample_number+1:n*sample_number);   
            en_1=imag(sum(r_T.*conj(ss))); 
            fi_est=fi_est+alpha*en_1;                          % 更新载波相位估计值 
            %@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
          end %else end 
          sum_Uk_old=sum_Uk_new;
          %******************************************************************************************************************* 
        end   % for n=1:len end

        max=1;
        for i=1:state_number
            if metric_state(i,n)>metric_state(max,n)
               max=i;
            end
        end
        final_msg_decision=msg_decision(max,:);
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        
        final_bit_decision=zeros(1,length(final_msg_decision)*2);
        for i=1:length(final_msg_decision)
            if final_msg_decision(i)==0
               final_bit_decision(2*i-1)=0;
               final_bit_decision(2*i)=0;
            elseif final_msg_decision(i)==1
               final_bit_decision(2*i-1)=0;
               final_bit_decision(2*i)=1; 
            elseif final_msg_decision(i)==2
               final_bit_decision(2*i-1)=1;
               final_bit_decision(2*i)=0;
            else
               final_bit_decision(2*i-1)=1;
               final_bit_decision(2*i)=1; 
            end
        end

        [TCM_bernum,TCM_berrat]=symerr(bit_msg,final_bit_decision);        % bit error rate
        BER(it)=BER(it-1)+TCM_bernum;
    end
    SNR
    BER_rat=BER(it)/it_num/40
    BER_SNR(tt)=BER_rat;
end
diary off

































